Tutorial

Go to the Download page to download this notebook

This is a basic tutorial on how to use ROSS (rotordynamics open-source software), a Python library for rotordynamic analysis. Most of this code follows object-oriented paradigm, which is represented in this  UML DIAGRAM.

Before starting the tutorial, it is worth noting some of ROSS’ design characteristics.

First, we can divide the use of ROSS in two steps: - Building the model; - Calculating the results.

We can build a model by instantiating elements such as beams (shaft), disks and bearings. These elements are all defined in classes with names such as ShaftElement, BearingElement and so on.

After instantiating some elements, we can then use these to build a rotor.

To get results, we always have to use one of the .run_ methods available for a rotor object. These methods will return objects that store the analysis results and that also have plot methods available. These methods will use the plotly library to make graphs common to a rotordynamic analysis.

Another important characteristic of ROSS, is that it uses the pint library to convert units.

This means that every time we call a function, we can use pint.Quantity as an argument for values that have units. If we give a float to the function ROSS will consider SI units as default.

We can also use units when plotting results. For example, for a unbalance response plot we have the amplitude_units argument and we can choose between any length unit available in pint such as ‘meter’, ‘inch’, etc.

In the following topics, we will discuss the most relevant classes for a quick start on how to use ROSS.

Table of Contents

• 1. Material Class

  • 1.1. Creating a material

  • 1.2. Saving materials

  • 1.3. Available materials

  • 1.4. Loading materials

• 2. ShaftElement Class

  • 2.1. Creating shaft elements

    • 2.1.1. Cylindrical shaft element

    • 2.1.2. Conical shaft element

    • 2.1.3. List of elements - identical properties

    • 2.1.4. List of elements - different properties

  • 2.2. Creating shaft elements via Excel

• 3. DiskElement Class

  • 3.1. Creating disk elements from inertia properties

    • 3.1.1. Creating a single disk element

    • 3.1.2. Creating a list of disk element

  • 3.2. Creating disk elements from geometrical properties

    • 3.2.1. Creating a single disk element

    • 3.2.2. Creating a list of disk element

  • 3.3. Creating disk elements via Excel

• 4. Bearing and Seal Classes

  • 4.1. BearingElement Class

    • 4.1.1. Bearing with constant coefficients

    • 4.1.2. Bearing with varying coefficients

    • 4.1.3. Inserting bearing elements in series

    • 4.1.4. Visualizing coefficients graphically

  • 4.2. SealElement Class

  • 4.3. BallBearingElement Class

  • 4.4. RollerBearingElement Class

  • 4.5. MagneticElement Class

    • 4.5.1. MagneticBearing from electromagnetic parameters

  • 4.6. Creating bearing elements via Excel

• 5. PointMass Class

• 6. Rotor Class

  • 6.1. Creating a rotor model

  • 6.2. Creating a rotor from sections

  • 6.3. Visualizing the rotor model

  • 6.4. Saving a rotor model

  • 6.5. Loading a rotor model

• 7. ROSS Units System

  • 7.1 Inserting units

  • 7.2 Manipulating units for plotting

• 8. Rotor Analyses

  • 8.1. Static Analysis

    • 8.1.1 .Running static analysis

    • 8.1.2. Plotting results

  • 8.2. Modal Analysis

    • 8.2.1. Running modal analysis

    • 8.2.2. Plotting results

  • 8.3. Campbell Diagram

    • 8.3.1. Running campbell diagram

    • 8.3.2. Plotting results

  • 8.4. Frequency Response

    • 8.4.1. Clustering points

    • 8.4.2. Running frequency response

    • 8.4.3. Plotting results - Bode Plot

    • 8.4.4. Plotting results - Deflected shape

  • 8.5. Unbalance Response

    • 8.5.1. Running unbalance response

    • 8.5.2. Plotting results - Bode Plot

    • 8.5.3. Plotting results - Deflected shape

  • 8.6. Time Response

    • 8.6.1. Running time response

    • 8.6.2. Plotting results

  • 8.7. Undamped Critical Speed Map (UCS)

    • 8.7.1. Running UCS

import os
from pathlib import Path
import ross as rs
import numpy as np

Section 1: Material Class

There is a class called Material to hold material’s properties. Materials are applied to shaft and disk elements.

1.1 Creating a material

To instantiate a Material class, you only need to give 2 out of the following parameters: E (Young’s Moduluds), G_s (Shear Modulus) ,Poisson (Poisson Coefficient), and the material density rho.

# from E and G_s
steel = rs.Material(name="Steel", rho=7810, E=211e9, G_s=81.2e9)
# from E and Poisson
steel2 = rs.Material(name="Steel", rho=7810, E=211e9, Poisson=0.3)
# from G_s and Poisson
steel3 = rs.Material(name="Steel", rho=7810, G_s=81.2e9, Poisson=0.3)

print(steel)

# returning attributes
print("="*36)
print(f"Young's Modulus: {steel.E}")
print(f"Shear Modulus:    {steel.G_s}")

Steel
-----------------------------------
Density         (kg/m**3): 7810.0
Young`s modulus (N/m**2):  2.11e+11
Shear modulus   (N/m**2):  8.12e+10
Poisson coefficient     :  0.29926108
====================================
Young's Modulus: 211000000000.0
Shear Modulus:    81200000000.0

Note: Adding 3 arguments to the Material class raises an error.

1.2 Saving materials

To save an already instantiated Material object, you need to use the following method.

steel.save_material()

1.3 Available materials

Saved Materials are stored in a .toml file, which can be read as .txt. The file is placed on ROSS root file with name available_materials.toml.

It’s possible to access the Material data from the file. With the file opened, you can: - modify the properties directly; - create new materials;

It’s important to keep the file structure to ensure the correct functioning of the class.

[Materials.Steel]
name = "Steel"
rho = 7810
E = 211000000000.0
Poisson = 0.2992610837438423
G_s = 81200000000.0
color = "#525252"

Do not change the dictionary keys and the order they’re built.

To check what materials are available, use the command:

rs.Material.available_materials()

['Steel', 'AISI4140', 'A216WCB']

1.4 Loading materials

After checking the available materials, you should use the Material.use_material('name') method with the name of the material as a parameter.

steel5 = rs.Material.load_material('Steel')

Section 2: ShaftElement Class

ShaftElement allows you to create cylindrical and conical shaft elements. It means you can set differents outer and inner diameters for each element node.

There are some ways in which you can choose the parameters to model this element:

• Euler–Bernoulli beam Theory (rotary_inertia=False, shear_effects=False)

• Timoshenko beam Theory (rotary_inertia=True, shear_effects=True - used as default)

The matrices (mass, stiffness, damping and gyroscopic) will be defined considering the following local coordinate vector:

\([x_0, y_0, \alpha_0, \beta_0, x_1, y_1, \alpha_1, \beta_1]^T\) Where \(\alpha_0\) and \(\alpha_1\) are the bending on the yz plane \(\beta_0\) and \(\beta_1\) are the bending on the xz plane.

This element represents the rotor’s shaft, all the other elements are correlated with this one.

2.1 Creating shaft elements

The next examples present different ways of how to create a ShaftElement object, from a single element to a list of several shaft elements with different properties.

When creating shaft elements, you don’t necessarily need to input a specific node. If n=None, the Rotor class will assign a value to the element when building a rotor model (see section 6).

You can also pass the same n value to several shaft elements in the same rotor model.

2.1.1 Cylindrical shaft element

As it’s been seen, a shaft element has 4 parameters for diameters. To simplify that, when creating a cylindrical element, you only need to give 2 of them: idl and odl. So the other 2 (idr and odr) get the same values.

Note: you can give all the 4 parameters, as long they match each other.

# Cylindrical shaft element
L = 0.25
i_d = 0
o_d = 0.05
cy_elem = rs.ShaftElement(
    L=L,
    idl=i_d,
    odl=o_d,
    material=steel,
    shear_effects=True,
    rotary_inertia=True,
    gyroscopic=True,
)

print(cy_elem)

Element Number:             None
Element Lenght   (m):       0.25
Left Int. Diam.  (m):        0.0
Left Out. Diam.  (m):       0.05
Right Int. Diam. (m):        0.0
Right Out. Diam. (m):       0.05
-----------------------------------
Steel
-----------------------------------
Density         (kg/m**3): 7810.0
Young`s modulus (N/m**2):  2.11e+11
Shear modulus   (N/m**2):  8.12e+10
Poisson coefficient     :  0.29926108

2.1.2 Conical shaft element

To create a conical shaft elements, you must give all the 4 diamater parameters, and idl != idr and/or odl != odr.

# Conical shaft element
L = 0.25
idl = 0
idr = 0
odl = 0.05
odr = 0.07
co_elem = rs.ShaftElement(
    L=L,
    idl=idl,
    idr=idr,
    odl=odl,
    odr=odr,
    material=steel,
    shear_effects=True,
    rotary_inertia=True,
    gyroscopic=True,
)
print(co_elem)

Element Number:             None
Element Lenght   (m):       0.25
Left Int. Diam.  (m):        0.0
Left Out. Diam.  (m):       0.05
Right Int. Diam. (m):        0.0
Right Out. Diam. (m):       0.07
-----------------------------------
Steel
-----------------------------------
Density         (kg/m**3): 7810.0
Young`s modulus (N/m**2):  2.11e+11
Shear modulus   (N/m**2):  8.12e+10
Poisson coefficient     :  0.29926108

Returning element matrices

Use one of this methods to return the matrices: - .M(): returns the mass matrix - .K(frequency): returns the stiffness matrix - .C(frequency): returns the damping matrix - .G(): returns de gyroscopic matrix

# Mass matrix
# cy_elem.M()

# Stiffness matrix
# frequency = 0
# cy_elem.K(frequency)

# Damping matrix
# frequency = 0
# cy_elem.C(frequency)

# Gyroscopic matrix
# cy_elem.G()

2.1.3 List of elements - identical properties

Now we leanrt how to create elements, let’s automate the process of creating multiple elements with identical properties.

In this example, we want 6 shaft elements with identical properties. This process can be done using a for loop or a list comprehension.

# Creating a list of shaft elements
L = 0.25
i_d = 0
o_d = 0.05
N = 6
shaft_elements = [
    rs.ShaftElement(
        L=L,
        idl=i_d,
        odl=o_d,
        material=steel,
        shear_effects=True,
        rotary_inertia=True,
        gyroscopic=True,
    )
    for _ in range(N)
]
shaft_elements

[ShaftElement(L=0.25, idl=0.0, idr=0.0, odl=0.05,  odr=0.05, material='Steel', n=None),
 ShaftElement(L=0.25, idl=0.0, idr=0.0, odl=0.05,  odr=0.05, material='Steel', n=None),
 ShaftElement(L=0.25, idl=0.0, idr=0.0, odl=0.05,  odr=0.05, material='Steel', n=None),
 ShaftElement(L=0.25, idl=0.0, idr=0.0, odl=0.05,  odr=0.05, material='Steel', n=None),
 ShaftElement(L=0.25, idl=0.0, idr=0.0, odl=0.05,  odr=0.05, material='Steel', n=None),
 ShaftElement(L=0.25, idl=0.0, idr=0.0, odl=0.05,  odr=0.05, material='Steel', n=None)]

2.1.4 List of elements - different properties

Now we leanrt how to create elements, let’s automate the process of creating multiple elements with identical properties.

In this example, we want 6 shaft elements which properties may not be the same. This process can be done using a for loop or a list comprehension, coupled with Python’s zip() method.

We create lists for each property, where each term refers to a single element:

# OPTION No.1:
# Using zip() method
L =   [0.20, 0.20, 0.10, 0.10, 0.20, 0.20]
i_d = [0.01,    0,    0,    0,    0, 0.01]
o_d = [0.05, 0.05, 0.06, 0.06, 0.05, 0.05]
shaft_elements = [
    rs.ShaftElement(
        L=l,
        idl=idl,
        odl=odl,
        material=steel,
        shear_effects=True,
        rotary_inertia=True,
        gyroscopic=True,
    )
    for l, idl, odl in zip(L, i_d, o_d)
]
shaft_elements

[ShaftElement(L=0.2, idl=0.01, idr=0.01, odl=0.05,  odr=0.05, material='Steel', n=None),
 ShaftElement(L=0.2, idl=0.0, idr=0.0, odl=0.05,  odr=0.05, material='Steel', n=None),
 ShaftElement(L=0.1, idl=0.0, idr=0.0, odl=0.06,  odr=0.06, material='Steel', n=None),
 ShaftElement(L=0.1, idl=0.0, idr=0.0, odl=0.06,  odr=0.06, material='Steel', n=None),
 ShaftElement(L=0.2, idl=0.0, idr=0.0, odl=0.05,  odr=0.05, material='Steel', n=None),
 ShaftElement(L=0.2, idl=0.01, idr=0.01, odl=0.05,  odr=0.05, material='Steel', n=None)]

# OPTION No.2:
# Using list index
L =   [0.20, 0.20, 0.10, 0.10, 0.20, 0.20]
i_d = [0.01,    0,    0,    0,    0, 0.01]
o_d = [0.05, 0.05, 0.06, 0.06, 0.05, 0.05]
N = len(L)
shaft_elements = [
    rs.ShaftElement(
        L=L[i],
        idl=i_d[i],
        odl=o_d[i],
        material=steel,
        shear_effects=True,
        rotary_inertia=True,
        gyroscopic=True,
    )
    for i in range(N)
]
shaft_elements

[ShaftElement(L=0.2, idl=0.01, idr=0.01, odl=0.05,  odr=0.05, material='Steel', n=None),
 ShaftElement(L=0.2, idl=0.0, idr=0.0, odl=0.05,  odr=0.05, material='Steel', n=None),
 ShaftElement(L=0.1, idl=0.0, idr=0.0, odl=0.06,  odr=0.06, material='Steel', n=None),
 ShaftElement(L=0.1, idl=0.0, idr=0.0, odl=0.06,  odr=0.06, material='Steel', n=None),
 ShaftElement(L=0.2, idl=0.0, idr=0.0, odl=0.05,  odr=0.05, material='Steel', n=None),
 ShaftElement(L=0.2, idl=0.01, idr=0.01, odl=0.05,  odr=0.05, material='Steel', n=None)]

2.2 Creating shaft elements via Excel

There is an option for creating a list of shaft elements via an Excel file. The classmethod .from_table() reads an Excel file created and converts it to a list of shaft elements.

A header with the names of the columns is required. These names should match the names expected by the routine (usually the names of the parameters, but also similar ones). The program will read every row bellow the header until they end or it reaches a NaN, which means if the code reaches to an empty line, it stops iterating.

An example of Excel content can be found at ROSS GitHub repository at ross/tests/data/shaft_si.xls, spreadsheet “Model”.

You can load it using the following code.

shaft_file = Path("shaft_si.xls")
shaft = rs.ShaftElement.from_table(
    file=shaft_file, sheet_type="Model", sheet_name="Model"
)

Section 3: DiskElement Class

The class DiskElement allows you to create disk elements, representing rotor equipments which can be considered only to add mass and inertia to the system, disregarding the stiffness.

ROSS offers 3 (three) ways to create a disk element: 1. Inputing mass and inertia data 2. Inputing geometrical and material data 3. From Excel table

3.1 Creating disk elements from inertia properties

If you have access to the mass and inertia properties of a equipment, you can input the data directly to the element.

Disk elements are useful to represent equipments which mass and inertia are significant, but the stiffness can be neglected.

3.1.1 Creating a single disk element

This example below shows how to instantiate a disk element according to the mass and inertia properties:

disk = rs.DiskElement(
    n=0,
    m=32.58,
    Ip=0.178,
    Id=0.329,
    tag="Disk"
)
disk

DiskElement(Id=0.329, Ip=0.178, m=32.58, color='Firebrick', n=0, tag='Disk')

3.1.2 Creating a list of disk element

This example below shows how to create a list of disk element according to the mass and inertia properties. The logic is the same applied to shaft elements.

# OPTION No.1:
# Using zip() method
n_list = [2, 4]
m_list = [32.6, 35.8]
Id_list = [0.17808928, 0.17808928]
Ip_list = [0.32956362, 0.38372842]
disk_elements = [
    rs.DiskElement(
        n=n,
        m=m,
        Id=Id,
        Ip=Ip,
    )
    for n, m, Id, Ip in zip(n_list, m_list, Id_list, Ip_list)
]
disk_elements

[DiskElement(Id=0.17809, Ip=0.32956, m=32.6, color='Firebrick', n=2, tag=None),
 DiskElement(Id=0.17809, Ip=0.38373, m=35.8, color='Firebrick', n=4, tag=None)]

# OPTION No.2:
# Using list index
n_list = [2, 4]
m_list = [32.6, 35.8]
Id_list = [0.17808928, 0.17808928]
Ip_list = [0.32956362, 0.38372842]
N = len(n_list)
disk_elements = [
    rs.DiskElement(
        n=n_list[i],
        m=m_list[i],
        Id=Id_list[i],
        Ip=Ip_list[i],
    )
    for i in range(N)
]
disk_elements

[DiskElement(Id=0.17809, Ip=0.32956, m=32.6, color='Firebrick', n=2, tag=None),
 DiskElement(Id=0.17809, Ip=0.38373, m=35.8, color='Firebrick', n=4, tag=None)]

3.2 Creating disk elements from geometrical properties

Besides the instantiation previously explained, there is a way to instantiate a DiskElement with only geometrical parameters (an approximation for cylindrical disks) and the disk’s material, as we can see in the following code. In this case, there’s a class method (rs.DiskElement.from_geometry()) which you can use.

ROSS will take geometrical parameters (outer and inner diameters, and width) and convert them into mass and inertia data. Once again, considering the disk as a cylinder.

3.2.1 Creating a single disk element

This example below shows how to instantiate a disk element according to the geometrical and material properties:

disk1 = rs.DiskElement.from_geometry(
    n=4,
    material=steel,
    width=0.07,
    i_d=0.05,
    o_d=0.28
)
print(disk1)

print("="*76)
print(f"Disk mass:              {disk1.m}")
print(f"Disk polar inertia:     {disk1.Ip}")
print(f"Disk diametral inertia: {disk1.Id}")

Tag:                      None
Node:                     4
Mass           (kg):      32.59
Diam. inertia  (kg*m**2): 0.17809
Polar. inertia (kg*m**2): 0.32956
============================================================================
Disk mass:              32.58972765304033
Disk polar inertia:     0.32956362089137037
Disk diametral inertia: 0.17808928257067666

3.2.2 Creating a list of disk element

This example below shows how to create a list of disk element according to the geometrical and material properties. The logic is the same applied to shaft elements.

# OPTION No.1:
# Using zip() method
n_list = [2, 4]
width_list = [0.7, 0.7]
i_d_list = [0.05, 0.05]
o_d_list = [0.15, 0.18]
disk_elements = [
    rs.DiskElement.from_geometry(
        n=n,
        material=steel,
        width=width,
        i_d=i_d,
        o_d=o_d,
    )
    for n, width, i_d, o_d in zip(n_list, width_list, i_d_list, o_d_list)
]
disk_elements

[DiskElement(Id=3.6408, Ip=0.26836, m=85.875, color='Firebrick', n=2, tag=None),
 DiskElement(Id=5.5224, Ip=0.56007, m=128.38, color='Firebrick', n=4, tag=None)]

# OPTION No.2:
# Using list index
n_list = [2, 4]
width_list = [0.7, 0.7]
i_d_list = [0.05, 0.05]
o_d_list = [0.15, 0.18]
N = len(n_list)
disk_elements = [
    rs.DiskElement.from_geometry(
        n=n_list[i],
        material=steel,
        width=width_list[i],
        i_d=i_d_list[i],
        o_d=o_d_list[i],
    )
    for i in range(N)
]
disk_elements

[DiskElement(Id=3.6408, Ip=0.26836, m=85.875, color='Firebrick', n=2, tag=None),
 DiskElement(Id=5.5224, Ip=0.56007, m=128.38, color='Firebrick', n=4, tag=None)]

3.3 Creating disk elements via Excel

The third option for creating disk elements is via an Excel file. The classmethod .from_table() reads an Excel file created and converts it to a list of disk elements. This method accepts only mass and inertia inputs.

A header with the names of the columns is required. These names should match the names expected by the routine (usually the names of the parameters, but also similar ones). The program will read every row bellow the header until they end or it reaches a NaN, which means if the code reaches to an empty line, it stops iterating.

You can take advantage of the excel file used to assemble shaft elements, to assemble disk elements, just add a new spreadsheet to your Excel file and specify the correct sheet_name.

An example of Excel content can be found at diretory ross/tests/data/shaft_si.xls, spreadsheet “More”.

file_path = Path("shaft_si.xls")
list_of_disks = rs.DiskElement.from_table(file=file_path, sheet_name="More")
list_of_disks

[DiskElement(Id=0.0, Ip=0.0, m=15.12, color='Firebrick', n=3, tag=None),
 DiskElement(Id=0.025, Ip=0.047, m=6.91, color='Firebrick', n=20, tag=None),
 DiskElement(Id=0.025, Ip=0.047, m=6.93, color='Firebrick', n=23, tag=None),
 DiskElement(Id=0.025, Ip=0.048, m=6.95, color='Firebrick', n=26, tag=None),
 DiskElement(Id=0.025, Ip=0.048, m=6.98, color='Firebrick', n=29, tag=None),
 DiskElement(Id=0.025, Ip=0.048, m=6.94, color='Firebrick', n=32, tag=None),
 DiskElement(Id=0.025, Ip=0.048, m=6.96, color='Firebrick', n=35, tag=None)]

Section 4: Bearing and Seal Classes

ROSS has a serie of classe to represent element that adds stiffness and / or damping to a rotor system. They’re suitable to represent mainly bearings, supports and seals. Each one aims to represent some types of bearing and seal.

All the class will return four stiffness coefficients (\(k_{xx}\), \(k_{xy}\), \(k_{yx}\), \(k_{yy}\)) and four damping coefficients (\(c_{xx}\), \(c_{xy}\), \(c_{yx}\), \(c_{yy}\)), which will be used to assemble the stiffness and damping matrices.

The main difference between these classes are the arguments the user must input to create the element.

Available bearing classes and class methods:

• 1. BearingElement: represents a general (journal) bearing element.

• 2. SealElement: represents a general seal element.

• 3. BallBearingElement: A bearing element for ball bearings

• 4. RollerBearingElement: A bearing element for roller bearings.

• 5. MagneticBearingElement: A bearing element for magnetic bearings.

  • 5.1. param_to_coef: A bearing element for magnetic bearings from electromagnetic parameters

The classes from item 2 to 5 inherits from BearingElement class. It means, you can use the same methods and commands, set up to BearingElement, in the other classes.

4.1 BearingElement Class

This class will create a bearing element. Bearings are elements that only add stiffness and damping properties to the rotor system. These parameters are defined by 8 dynamics coefficients (4 stiffness coefficients and 4 damping coefficients).

Parameters can be a constant value or speed dependent. For speed dependent parameters, each argument should be passed as an array and the correspondent speed values should also be passed as an array. Values for each parameter will be interpolated for the speed.

Bearing elements are single node elements and linked to “ground”, but it’s possible to create a new node with n_link argument to introduce a link with other elements. Useful to add bearings in series or co-axial rotors.

4.1.1 Bearing with constant coefficients

Bearings can have a constant value for each coefficient. In this case, it’s not necessary to give a value to frequency argument.

The next example shows how to instantiate a single bearing with constant coefficients:

stfx = 1e6
stfy = 0.8e6
bearing1 = rs.BearingElement(n=0, kxx=stfx, kyy=stfy, cxx=1e3)
print(bearing1)

print("="*55)
print(f"Kxx coefficient - array of strings format: {bearing1.kxx}")
print(f"Kxx coefficient - array of floats format:  {bearing1.kxx.coefficient}")

BearingElement(n=0, n_link=None,
 kxx=[1000000.0], kxy=[0],
 kyx=[0], kyy=[800000.0],
 cxx=[1000.0], cxy=[0],
 cyx=[0], cyy=[1000.0],
 frequency=None, tag=None)
=======================================================
Kxx coefficient - array of strings format: [1000000.0]
Kxx coefficient - array of floats format:  [1000000.0]

4.1.2 Bearing with varying coefficients

The coefficients could be an array with different values for different rotation speeds, in that case you only have to give a parameter ‘frequency’ which is a array with the same size as the coefficients array.

The next example shows how to instantiate a single bearing with speed dependent parameters:

bearing2 = rs.BearingElement(
    n=0,
    kxx=np.array([0.5e6, 1.0e6, 2.5e6]),
    kyy=np.array([1.5e6, 2.0e6, 3.5e6]),
    cxx=np.array([0.5e3, 1.0e3, 1.5e3]),
    frequency=np.array([0, 1000, 2000]),
)

print(bearing2)
print("="*79)
print(f"Kxx coefficient - array of strings format: {bearing2.kxx}")
print(f"Kxx coefficient - array of floats format:  {bearing2.kxx.coefficient}")

BearingElement(n=0, n_link=None,
 kxx=array([ 500000., 1000000., 2500000.]), kxy=[0, 0, 0],
 kyx=[0, 0, 0], kyy=array([1500000., 2000000., 3500000.]),
 cxx=array([ 500., 1000., 1500.]), cxy=[0, 0, 0],
 cyx=[0, 0, 0], cyy=array([ 500., 1000., 1500.]),
 frequency=[   0. 1000. 2000.], tag=None)
===============================================================================
Kxx coefficient - array of strings format: array([ 500000., 1000000., 2500000.])
Kxx coefficient - array of floats format:  [ 500000. 1000000. 2500000.]

If the size of coefficient and frequency arrays do not match, an ValueError is raised

The next example shows the instantiate of a bearing with odd parameters:

bearing_odd = rs.BearingElement( # odd dimensions
    n=0,
    kxx=np.array([0.5e6, 1.0e6, 2.5e6]),
    kyy=np.array([1.5e6, 2.0e6, 3.5e6]),
    cxx=np.array([0.5e3, 1.0e3, 1.5e3]),
    frequency=np.array([0, 1000, 2000, 3000])
)

4.1.3 Inserting bearing elements in series

Bearing and seal elements are 1-node element, which means the element attaches to a given node from the rotor shaft and it’s connect to the “ground”. However, there’s an option to couple multiple elements in series, using the n_link argument. This is very useful to simulate structures which support the machine, for example.

n_link opens a new node to the rotor system, or it can be associated to another rotor node (useful in co-axial rotor models). Then, the new BearingElement node, is set equal to the n_link from the previous element.

stfx = 1e6
stfy = 0.8e6
bearing3 = rs.BearingElement(n=0, kxx=stfx, kyy=stfy, cxx=1e3, n_link=1, tag="journal_bearing")
bearing4 = rs.BearingElement(n=1, kxx=1e7, kyy=1e9, cxx=10, tag="support")
print(bearing3)
print(bearing4)

BearingElement(n=0, n_link=1,
 kxx=[1000000.0], kxy=[0],
 kyx=[0], kyy=[800000.0],
 cxx=[1000.0], cxy=[0],
 cyx=[0], cyy=[1000.0],
 frequency=None, tag='journal_bearing')
BearingElement(n=1, n_link=None,
 kxx=[10000000.0], kxy=[0],
 kyx=[0], kyy=[1000000000.0],
 cxx=[10], cxy=[0],
 cyx=[0], cyy=[10],
 frequency=None, tag='support')

4.1.4 Visualizing coefficients graphically

If you want to visualize how the coefficients varies with speed, you can select a specific coefficient and use the .plot() method.

Let’s return to the example done in 4.1.2 and check how \(k_{yy}\) and \(c_{yy}\) varies. You can check for all the 8 dynamic coefficients as you like.

bearing2.kyy.plot()

bearing2.cyy.plot()

4.2 SealElement Class

SealElement class method have the exactly same arguments than BearingElement. The differences are found in some considerations when assembbling a full rotor model. For example, a SealElement won’t generate reaction forces in a static analysis. So, even they are very similar when built, they have different roles in the model.

Let’s see an example:

stfx = 1e6
stfy = 0.8e6
seal = rs.SealElement(n=0, kxx=stfx, kyy=stfy, cxx=1e3, cyy=0.8e3)
seal

SealElement(n=0, n_link=None,
 kxx=[1000000.0], kxy=[0],
 kyx=[0], kyy=[800000.0],
 cxx=[1000.0], cxy=[0],
 cyx=[0], cyy=[800.0],
 frequency=None, tag=None)

4.3 BallBearingElement Class

This class will create a bearing element based on some geometric and constructive parameters of ball bearings. The main difference is that cross-coupling stiffness and damping are not modeled in this case.

Let’s see an example:

n = 0
n_balls= 8
d_balls = 0.03
fs = 500.0
alpha = np.pi / 6
tag = "ballbearing"
ballbearing = rs.BallBearingElement(
    n=n,
    n_balls=n_balls,
    d_balls=d_balls,
    fs=fs,
    alpha=alpha,
    tag=tag,
)
ballbearing

BallBearingElement(n=0, n_link=None,
 kxx=[46416883.847697675], kxy=[0.0],
 kyx=[0.0], kyy=[100906269.23412538],
 cxx=[580.211048096221], cxy=[0.0],
 cyx=[0.0], cyy=[1261.3283654265672],
 frequency=None, tag='ballbearing')

4.4 RollerBearingElement Class

This class will create a bearing element based on some geometric and constructive parameters of roller bearings. The main difference is that cross-coupling stiffness and damping are not modeled in this case.

Let’s see an example:

n = 0
n_rollers= 8
l_rollers = 0.03
fs = 500.0
alpha = np.pi / 6
tag = "rollerbearing"
rollerbearing = rs.RollerBearingElement(
    n=n,
    n_rollers=n_rollers,
    l_rollers=l_rollers,
    fs=fs,
    alpha=alpha,
    tag=tag
)
rollerbearing

RollerBearingElement(n=0, n_link=None,
 kxx=[272821927.4006065], kxy=[0.0],
 kyx=[0.0], kyy=[556779443.6747072],
 cxx=[3410.2740925075814], cxy=[0.0],
 cyx=[0.0], cyy=[6959.74304593384],
 frequency=None, tag='rollerbearing')

4.5 MagneticBearingElement Class

This class creates a magnetic bearing element. You can input electromagnetic parameters and PID gains. ROSS converts it to stiffness and damping coefficients. To do ir, use the class MagneticBearingElement()

See the following reference for the electromagnetic parameters g0, i0, ag, nw, alpha: Book: Magnetic Bearings. Theory, Design, and Application to Rotating Machinery Authors: Gerhard Schweitzer and Eric H. Maslen Page: 84-95

From: “Magnetic Bearings. Theory, Design, and Application to Rotating Machinery” Authors: Gerhard Schweitzer and Eric H. Maslen Page: 354

Let’s see an example:

n = 0
g0 = 1e-3
i0 = 1.0
ag = 1e-4
nw = 200
alpha = 0.392
kp_pid = 1.0
kd_pid = 1.0
k_amp = 1.0
k_sense = 1.0
tag = "magneticbearing"
mbearing = rs.MagneticBearingElement(
    n=n,g0=g0,i0=i0,ag=ag,nw=nw,alpha=alpha, kp_pid=kp_pid,kd_pid=kd_pid, k_amp=k_amp, k_sense=k_sense, tag=tag
)
mbearing

MagneticBearingElement(n=0, n_link=None,
 kxx=[-4640.623377181318], kxy=[0.0],
 kyx=[0.0], kyy=[-4640.623377181318],
 cxx=[4.645268645827145], cxy=[0.0],
 cyx=[0.0], cyy=[4.645268645827145],
 frequency=None, tag='magneticbearing')

4.6 Creating bearing elements via Excel

There’s an option for creating bearing elements via an Excel file. The classmethod .from_table() reads an Excel file created and converts it to a BearingElement instance. Differently from creating shaft or disk elements, this method creates only a single bearing element. To create a list of bearing elements, the user should open several spreadsheets in the Excel file and run a list comprehension loop appending each elemnet to the list.

A header with the names of the columns is required. These names should match the names expected by the routine (usually the names of the parameters, but also similar ones). The program will read every row bellow the header until they end or it reaches a NaN, which means if the code reaches to an empty line, it stops iterating.

n : int
    The node in which the bearing will be located in the rotor.
file: str
    Path to the file containing the bearing parameters.
sheet_name: int or str, optional
    Position of the sheet in the file (starting from 0) or its name. If none is passed, it is
    assumed to be the first sheet in the file.

An example of Excel content can be found at diretory ross/tests/data/bearing_seal_si.xls, spreadsheet “XLUserKCM”.

# single bearing element
file_path = Path("bearing_seal_si.xls")
bearing = rs.BearingElement.from_table(n=0, file=file_path)
bearing

BearingElement(n=0, n_link=None,
 kxx=array([1.37981e+07, 2.99519e+07, 5.35657e+07, 8.51442e+07, 1.20733e+08,
       1.59519e+08, 1.97885e+08, 2.35240e+08, 2.71250e+08]), kxy=array([0, 0, 0, 0, 0, 0, 0, 0, 0]),
 kyx=array([0, 0, 0, 0, 0, 0, 0, 0, 0]), kyy=array([1.37981e+07, 2.99519e+07, 5.35657e+07, 8.51442e+07, 1.20733e+08,
       1.59519e+08, 1.97885e+08, 2.35240e+08, 2.71250e+08]),
 cxx=array([102506, 127450, 144989, 153563, 155122, 150835, 145086, 141871,
       140702]), cxy=array([0, 0, 0, 0, 0, 0, 0, 0, 0]),
 cyx=array([0, 0, 0, 0, 0, 0, 0, 0, 0]), cyy=array([102506, 127450, 144989, 153563, 155122, 150835, 145086, 141871,
       140702]),
 frequency=[ 314.15926536  418.87902048  523.5987756   628.31853072  733.03828584
  837.75804096  942.47779608 1047.1975512  1151.91730632], tag=None)

As .from_table() creates only a single bearing, let’s see an example how to create multiple elements without typing the same command line multiple times.

• First, in the EXCEL file, create multiple spreadsheets. Each one must hold the bearing coefficients and frequency data.

• Then, create a list holding the node numbers for each bearing (respecting the order of the spreadsheets from the EXCEL file).

• Finally, create a loop which iterates over the the nodes list and the spreadsheet.

# list of bearing elements

# nodes = list with the bearing elements nodes number
# file_path = Path("bearing_seal_si.xls")
# bearings = [rs.BearingElement.from_table(n, file_path, sheet_name=i) for i, n in enumerate(nodes)]

Section 5: PointMass Class

The PointMass class creates a point mass element. This element can be used to link other elements in the analysis. The mass provided can be different on the x and y direction (e.g. different support inertia for x and y directions).

PointMass also keeps the mass, stiffness, damping and gyroscopic matrices sizes consistence. When adding 2 bearing elements in series, it opens a new node with new degrees of freedom (DoF) (see section 4.1.3) and expands the stiffness and damping matrices. For this reason, it’s necessary to add mass values to those DoF to match the matrices sizes.

If you input the argument m, the code automatically replicate the mass value for both directions “x” and “y”.

Let’s see an example of creating point masses:

# inputting m
p0 = rs.PointMass(n=0, m=2)
p0.M() # returns de mass matrices for the element

array([[2., 0.],
       [0., 2.]])

# inputting mx and my
p1 = rs.PointMass(n=0, mx=2, my=3)
p1.M()

array([[2., 0.],
       [0., 3.]])

Section 6: Rotor Class

Rotor is the main class from ROSS. It takes as argument lists with all elements and assembles the mass, gyroscopic, damping and stiffness global matrices for the system. The object created has several methods that can be used to evaluate the dynamics of the model (they all start with the prefix .run_).

To use this class, you must input all the already instantiated elements in a list format.

If the shaft elements are not numbered, the class set a number for each one, according to the element’s position in the list supplied to the rotor constructor.

To assemble the matrices, the Rotor class takes the local DoF’s index from each element (element method .dof_mapping()) and calculate the global index

6.1 Creating a rotor model

Let’s create a simple rotor model with \(1.5 m\) length with 6 identical shaft elements, 2 disks, 2 bearings in the shaft ends and a support linked to the first bearing. First, we create the elements, then we input them to the Rotor class.

n = 6

shaft_elem = [
    rs.ShaftElement(
        L=0.25,
        idl=0.0,
        odl=0.05,
        material=steel,
        shear_effects=True,
        rotary_inertia=True,
        gyroscopic=True,
    )
    for _ in range(n)
]

disk0 = rs.DiskElement.from_geometry(
    n=2, material=steel, width=0.07, i_d=0.05, o_d=0.28
)
disk1 = rs.DiskElement.from_geometry(
    n=4, material=steel, width=0.07, i_d=0.05, o_d=0.28
)
disks = [disk0, disk1]

stfx = 1e6
stfy = 0.8e6
bearing0 = rs.BearingElement(0, kxx=stfx, kyy=stfy, cxx=0, n_link=7)
bearing1 = rs.BearingElement(6, kxx=stfx, kyy=stfy, cxx=0)
bearing2 = rs.BearingElement(7, kxx=stfx, kyy=stfy, cxx=0)

bearings = [bearing0, bearing1, bearing2]

pm0 = rs.PointMass(n=7, m=30)
pointmass = [pm0]

rotor1 = rs.Rotor(shaft_elem, disks, bearings, pointmass)

print("Rotor total mass = ", np.round(rotor1.m, 2))
print("Rotor center of gravity =", np.round(rotor1.CG, 2))

# plotting the rotor model
rotor1.plot_rotor()

Rotor total mass =  88.18
Rotor center of gravity = 0.75

6.2 Creating a rotor from sections

An alternative to build rotor models is dividing the rotor in sections. Each section gets the same number of shaft elements.

There’s an important difference in this class method when placing disks and bearings. The argument n will refer, not to the element node, but to the section node. So if your model has 3 sections with 4 elements each, there’re 4 section nodes and 13 element nodes.

Let’s repeat the rotor model from the last example, but using .from_section() class method, without the support.

i_d = 0
o_d = 0.05

# inner diameter of each section
i_ds_data = [0, 0, 0]
# outer diameter of each section
o_ds_data = [0.05, 0.05, 0.05]
# length of each section
leng_data = [0.5, 0.5, 0.5]

material_data = [steel, steel, steel]
# material_data = steel
stfx = 1e6
stfy = 0.8e6

# n = 0 refers to the section 0, first node
bearing0 = rs.BearingElement(n=0, kxx=stfx, kyy=stfy, cxx=1e3)

# n = 3 refers to the section 2, last node
bearing1 = rs.BearingElement(n=3, kxx=stfx, kyy=stfy, cxx=1e3)
bearings = [bearing0, bearing1]

# n = 1 refers to the section 1, first node
disk0 = rs.DiskElement.from_geometry(
    n=1,
    material=steel,
    width=0.07,
    i_d=0.05,
    o_d=0.28
)

# n = 2 refers to the section 2, first node
disk1 = rs.DiskElement.from_geometry(
    n=2,
    material=steel,
    width=0.07,
    i_d=0.05,
    o_d=0.28
)
disks = [disk0,disk1]

rotor2 = rs.Rotor.from_section(
    brg_seal_data=bearings,
    disk_data=disks,
    idl_data=i_ds_data,
    leng_data=leng_data,
    odl_data=o_ds_data,
    nel_r=4,
    material_data=steel,
)

print("Rotor total mass = ", np.round(rotor2.m, 2))
print("Rotor center of gravity =", np.round(rotor2.CG, 2))
rotor2.plot_rotor()

Rotor total mass =  88.18
Rotor center of gravity = 0.75

6.3 Visualizing the rotor model

It is interesting to plot the rotor to check if the geometry checks with what you wanted to the model. Use the .plot_rotor() method to create a plot.

nodes argument is useful when your model has lots of nodes and the visualization of nodes label may be confusing. Set an increment to the plot nodes label

ROSS uses PLOTLY as main plotting library:

With the Plotly, you can hover the mouse icon over the shaft, disk and point mass elements to check some of their parameters.

rotor1.plot_rotor()

Let’s visualize another rotor example with overlapping shaft elements:

shaft_file = Path("shaft_si.xls")
shaft = rs.ShaftElement.from_table(
    file=shaft_file, sheet_type="Model", sheet_name="Model"
)

file_path = Path("shaft_si.xls")
list_of_disks = rs.DiskElement.from_table(file=file_path, sheet_name="More")

bearing1 = rs.BearingElement.from_table(n=7, file="bearing_seal_si.xls", scale_factor=0.5)
bearing2 = rs.BearingElement.from_table(n=48, file="bearing_seal_si.xls", scale_factor=0.5)

bearings = [bearing1, bearing2]

rotor3 = rs.Rotor(shaft, list_of_disks, bearings)

node_increment = 5
rotor3.plot_rotor(nodes=node_increment)

6.4 Saving a rotor model

You can save a rotor model using the method .save(). This method saves the each element type and the rotor object in different .toml files.

You just need to input a name and the diretory, where it will be saved. If you don’t input a file_path, the rotor model is saved inside the “ross” folder.

To save the rotor2 we can use:

rotor2.save('rotor2.toml')

6.5 Loading a rotor model

You can load a rotor model using the method .load(). This method loads a previously saved rotor model.

You just need to input the file path to the method.

Now, let’s load the rotor2 we saved before:

rotor2_1 = rs.Rotor.load('rotor2.toml')
rotor2_1 == rotor2

True

Section 7: ROSS Units System

ROSS uses an units system package called Pint.

Pint defines, operates and manipulates physical quantities: the product of a numerical value and a unit of measurement. It allows arithmetic operations between them and conversions from and to different units.

With Pint, it’s possible to define units to every element type available in ROSS and manipulate the units when plotting graphs. ROSS takes the user-defined units and internally converts them to the International System (SI).

Important: It’s not possible to manipulate units for attributes from any class. Attributes’ values are always returned converted to SI. Only plot methods are able to manipulate the output unit.

7.1 Inserting units

Working with Pint requires a specific syntax to assign an unit to an argument.

First of all, it’s necessary to import a function called Q_ from ross.units. This function must be assigned to every variable that are desired to have units, followed by a tuple containing the magnitude and the unit (in string format).

The example below shows how to create a material using Pint, and how it is returned to the user.

from ross.units import Q_

rho = Q_(487.56237, "lb/foot**3")  # Imperial System
E = Q_(211.e9, "N/m**2")           # International System
G_s=Q_(81.2e9, "N/m**2")           # International System

steel4 = rs.Material(name="steel", rho=rho, E=E, G_s=G_s)

Note: Taking a closer look to the output values, the material density is converted to the SI and it’s returned this way to the user.

print(steel4)

steel
-----------------------------------
Density         (kg/m**3): 7810.0
Young`s modulus (N/m**2):  2.11e+11
Shear modulus   (N/m**2):  8.12e+10
Poisson coefficient     :  0.29926108

The same syntax applies to elements instantiation, if units are desired. Besides, notice the output is displayed in SI units.

Shaft Element using Pint

L = Q_(10, "in")
i_d = Q_(0., "meter")
o_d = Q_(0.05, "meter")

elem_pint = rs.ShaftElement(L=L, idl=i_d, odl=o_d, material=steel)
print(elem_pint)

Element Number:             None
Element Lenght   (m):      0.254
Left Int. Diam.  (m):        0.0
Left Out. Diam.  (m):       0.05
Right Int. Diam. (m):        0.0
Right Out. Diam. (m):       0.05
-----------------------------------
Steel
-----------------------------------
Density         (kg/m**3): 7810.0
Young`s modulus (N/m**2):  2.11e+11
Shear modulus   (N/m**2):  8.12e+10
Poisson coefficient     :  0.29926108

Bearing Element using Pint

kxx = Q_(2.54e4, "N/in")
cxx = Q_(1e2, "N*s/m")

brg_pint = rs.BearingElement(n=0, kxx=kxx, cxx=cxx)
print(brg_pint)

BearingElement(n=0, n_link=None,
 kxx=[1000000.0], kxy=[0],
 kyx=[0], kyy=[1000000.0],
 cxx=[100.0], cxy=[0],
 cyx=[0], cyy=[100.0],
 frequency=None, tag=None)

7.2 Manipulating units for plotting

The plot methods presents arguments to change the units for each axis. This kind of manipulation does not affect the resulting data stored. It only converts the data on the graphs.

The arguments names follow a simple logic. It is the “axis name” underscore “units” (axisname_units). It should help the user to identify which axis to modify. For example:

• frequency_units:

  • “rad/s”, “RPM”, “Hz”…

• amplitude_units:

  • “m”, “mm”, “in”, “foot”…

• displacement_units:

  • “m”, “mm”, “in”, “foot”…

• rotor_length_units:

  • “m”, “mm”, “in”, “foot”…

• moment_units:

  • “N/m”, “lbf/foot”…

It’s not necessary to add units previously to each element or material to use Pint with plots. But keep in mind ROSS will considers results values in the SI units.

Note: If you input data using the Imperial System, for example, without using Pint, ROSS will consider it’s in SI if you try to manipulate the units when plotting.

Let’s run a simple example of manipulating units for plotting.

samples = 31
speed_range = np.linspace(315, 1150, samples)

campbell = rotor3.run_campbell(speed_range)

Plotting with default options will bring graphs with SI units. X and Y axes representing the frequencies are set to rad/s

campbell.plot()

Now, let’s change the units to RPM.

Just by adding frequency_units="rpm" to plot method, you’ll change the plot units.

campbell.plot(frequency_units="RPM")

Section 8: Rotor Analyses

In section 6 we have learnt how to create a rotor model with Rotor class. Now, we’ll use the same class to run the simulation. There’re some methods, most of them with the prefix run_ you can use to run the rotordynamics analyses.

For Most of the methods, ypu can use the command .plot() to display a graphical visualization of the results (e.g run_campbel().plot(), run_freq_response().plot()).

ROSS offers the following analyses: - Static analysis - Modal analysis - Campbell Diagram - Frequency response - Unbalance response - Time response

Plotly library

ROSS uses Plotly for plotting results. All the figures can be stored and manipulated following Plotly API.

The following sections presents the results and how to return the Plotly Figures.

8.1 Static Analysis

This method runs the static analysis for the rotor. It calculate the static deformation due the gravity effects (shaft and disks weight). It also returns the bending moment and shearing force on each node, and you can return a free-body-diagram representation for the rotor, with the self weight, disks weight and reaction forces on bearings displayed.

8.1.1 Running static analysis

To run the simulation, use the .run_static() method. You can define a variable to store the results.

Storing the results, it’s possible to return the following arrays: - disk_forces_nodal - disk_forces_tag - bearing_forces_nodal - bearing_forces_tag - disp_y - Vx - Bm

static = rotor3.run_static()

Returning forces

Disk forces

• .disk_forces_nodal: Returns a dictionaty expliciting the node where the disk is located and the force value.

• .disk_forces_tag: Returns a dictionaty expliciting the the disk tag is located and the force value.

Bearing forces

• .bearing_forces_nodal: Returns a dictionaty expliciting the node where the bearing is located and the force value.

• .bearing_forces_tag: Returns a dictionaty expliciting the the bearing tag is located and the force value.

print("Disk forces - nodes")
print(rotor3.disk_forces_nodal)
print("")
print("Disk forces - tags")
print(rotor3.disk_forces_tag)

print("")
print("Bearing forces - nodes")
print(rotor3.bearing_forces_nodal)
print("")
print("Bearing forces - tags")
print(rotor3.bearing_forces_tag)

Disk forces - nodes
{'node_3': -148.3, 'node_20': -67.8, 'node_23': -68.0, 'node_26': -68.2, 'node_29': -68.4, 'node_32': -68.1, 'node_35': -68.3}

Disk forces - tags
{'Disk 0': -148.3, 'Disk 1': -67.8, 'Disk 2': -68.0, 'Disk 3': -68.2, 'Disk 4': -68.4, 'Disk 5': -68.1, 'Disk 6': -68.3}

Bearing forces - nodes
{'node_7': 1216.3, 'node_48': 1204.7}

Bearing forces - tags
{'Bearing 0': 1216.3, 'Bearing 1': 1204.7}

Other attributes from static analysis

• .Vx: Shearing force array

• .Bm: Bending moment array

• .deformation Displacement in Y direction

8.1.2 Plotting results

With results stored, you can use some methods to plot the results. Currently, there’re four plots you can retrieve from static analysis:

• .plot_free_body_diagram()

• .plot_deformation()

• .plot_shearing_force()

• .plot_bending_moment()

Plotting free-body-diagram

static.plot_free_body_diagram()

Plotting deformation

static.plot_deformation()

Plotting shearing force diagram

static.plot_shearing_force()

Plotting bending moment diagram

static.plot_bending_moment()

8.2 Modal Analysis

ROSS performs the modal analysis through method run_modal(). This method calculates natural frequencies, damping ratios and mode shapes.

You must select a speed, which will be used as excitation frequency to calculate the system’s eigenvalues and eigenvectors, and the number of eigenvalues and eigenvectors to be calculated is an optional argument (num_modes).

After running the modal analysis, it’s possible to return the following attributes: - eigenvalues (evalues); - eigenvectors (evectors); - damped natural frequencies (wd); - undamped natural frequencies (wn); - damping ratio (damping_ratio); - logarithmic decrement (log_dec).

8.2.1 Running modal analysis

To run the modal analysis, choose a speed to instantiate the method. For different speeds, change the the argument and run run_modal() once again.

Returning undamped natural frequencies

rotor_speed = 100.0 # rad/s
modal = rotor3.run_modal(rotor_speed)
print(f"Undamped natural frequencies:\n {modal.wn}")

Undamped natural frequencies:
 [ 87.53726224 470.55037053  96.28151251 724.61561509 910.9075301
 916.72762808]

Returning damped natural frequencies

# modal.wd
print(f"Damped natural frequencies:\n {modal.wd}")

Damped natural frequencies:
 [6.63919079e-05 1.60822814e-01 2.40535187e-01 1.77892662e+00
 9.02957212e+02 9.08665889e+02]

Returning the damping ratio

# modal.damping_ratio
print(f"Damping ratio for each mode:\n {modal.damping_ratio}")

Damping ratio for each mode:
 [-1.          0.99999994 -0.99999688  0.99999699  0.13183187  0.13232817]

Returning logarithmic decrement

# modal.log_dec
print(f"Logarithmic decrement for each mode:\n {modal.log_dec}")

Logarithmic decrement for each mode:
 [-8.28373218e+06  1.83839278e+04 -2.51502786e+03  2.55934136e+03
  8.35617240e-01  8.38819001e-01]

8.2.2 Plotting results

Once run_modal() is completed, you can check for the rotor’s mode shapes. You can plot each one of the modes calculated.

Besides, there’re two options for visualization: - plot_mode_2d - plotting 2D view - plot_mode_3d - plotting 3D view

Plotting 3D view

Use the command .plot_mode_3d(mode).

mode = 5
modal.plot_mode_3d(mode)

Plotting 2D view

Use the command .plot_mode_2d(mode).

mode = 5
modal.plot_mode_2d(mode)

8.3 Campbell Diagram

Also called “whirl speed map” in rotordynamics, ROSS calculate and plots the modes’ damped eigenvalues and the logarithmic decrement as a function of rotor speed.

To run the Campbell Diagram, use the command .run_campbell(). The user must input an array of speeds, which will be iterated to calculate each point on the graph.

This method returns the damped natural frequencies, logarithmic decrement and the whirl values (values indicating the whirl direction: backward or forward).

8.3.1 Running campbell diagram

In this example the whirl speed map is calculated for a speed range from 0 to 1000 rad/s (~9550 RPM).

Storing the results, it’s possible to return the following arrays: - wd - log_dec - whirl_values

Each value in these arrays is calculated for each speed value in speed_range

samples = 31
speed_range = np.linspace(315, 1150, samples)

campbell = rotor3.run_campbell(speed_range)

# results for each frequency
frequency_index = 0
print(campbell.wd[:, frequency_index])

[4.05696151e-03 1.13308765e-02 2.87443806e-02 6.81662768e-02
 1.53837839e-01 3.34102887e-01 7.04288351e-01 1.45788447e+00
 3.03663624e+00 5.25443046e+00 7.21720677e+00 1.06507330e+01
 1.85851449e+01 1.01602531e+02 4.67359814e+02 6.34570746e+02
 6.35176758e+02 6.35970906e+02 6.36891161e+02 6.37880817e+02
 6.38897245e+02 6.39918527e+02 6.40931104e+02 6.41925672e+02
 6.42889625e+02 6.43808455e+02 6.44672982e+02 6.45478063e+02
 6.46221340e+02 6.46902349e+02 6.47521893e+02]

# results for each frequency
frequency_index = 1
print(campbell.log_dec[:, frequency_index])

[1.29458650e+03 1.24034622e+03 1.18625615e+03 1.13071169e+03
 1.07242137e+03 1.01021334e+03 9.42865194e+02 8.68904404e+02
 7.86381133e+02 7.82375743e+02 3.44920711e+02 3.94207681e+02
 2.88161979e+02 7.77554137e+01 1.82932326e+01 4.83302025e-01
 4.33157440e-01 3.88197688e-01 3.48380244e-01 3.13495797e-01
 2.83129231e-01 2.56672664e-01 2.33605438e-01 2.13500543e-01
 1.95967979e-01 1.80635634e-01 1.67181122e-01 1.55328894e-01
 1.44844484e-01 1.35528792e-01 1.27212823e-01]

8.3.2 Plotting results

Now that the results are stored, use .plot() method to display the Campbell Diagram plot.

For the Campbell Diagram, you can plot more than one harmonic. As default, the plot display only the 1x speed option. Input a list with more harmonics to display it at the graph.

campbell.plot(harmonics=[0.5, 1])

8.4 Frequency Response

ROSS’ method to calculate the Frequency Response Function is run_freq_response(). This method returns the magnitude and phase in the frequency domain. The response is calculated for each node from the rotor model.

When plotting the results, you can choose to plot: - amplitude vs frequency: plot_magnitude() - phase vs frequency: plot_phase() - polar plot of amplitude vs phase: plot_polar_bode() - all: plot()

8.4.1 Clustering points

The number of solution points is an important parameter to determine the computational cost of the simulation. Besides the classical method, using numpy.linspace, which creates an evenly spaced array over a specified interval, ROSS offers an automatic method to create an speed_range array.

The method clustering_points generates an automatic array to run frequency response analyses. The frequency points are calculated based on the damped natural frequencies and their respective damping ratios. The greater the damping ratio, the more spread the points are. If the damping ratio, for a given critical speed, is smaller than 0.005, it is redefined to be 0.005 (for this method only).

The main goal of this feature is getting a more accurate amplitude value for the respective critical frequencies and nearby frequency points.

8.4.2 Running frequency response

To run the this analysis, use the command run_freq_response(). You can give a specific speed_range or let the program run with the default options. In this case, no arguments are needed to input.

First, let’s run an example with a “user-defined” speed_range. Setting an array to speed_range will disable all the frequency spacing parameters.

samples = 61
speed_range = np.linspace(315, 1150, samples) # rads/s
results1 = rotor3.run_freq_response(speed_range=speed_range)

results1.speed_range.size

61

Now, let’s run an example using clustering points array.

# results1_2 = rotor2.run_freq_response(cluster_points=True, num_points=5, num_modes=12)

# results1_2.speed_range.size

In the next section we’ll check the difference between both results.

8.4.3 Plotting results - Bode Plot

We can plot the frequency response of selecting the input and output degree of freedom.

• Input is the degree of freedom to be excited;

• Output is the degree of freedom to be observed.

Each shaft node has 4 local degrees of freedom (dof) \([x, y, \alpha, \beta]\), and each degree of freedom has it own index: - \(x\) -> index 0 - \(y\) -> index 1 - \(\alpha\) -> index 2 - \(\beta\) -> index 3

To select a DoF to input and a DoF to the output, we have to use the following correlation:

\(global\_dof = dof\_per\_node \cdot node\_number + dof\_index\)

For example: node 26, local dof \(y\):

\(DoF = 26 \cdot 4 + 1 = 25\)

plot = results1.plot(inp=105, out=105)
plot

# converting the first plot yaxis to log scale
# plot = results1.plot(inp=105, out=105)
# plot.update_yaxes(type="log", row=1, col=1)
# plot

# plot1_2 = results1_2.plot(inp=105, out=105)
# plot1_2

# # converting the first plot yaxis to log scale
# plot1_2 = results1_2.plot(inp=105, out=105)
# plot1_2.update_yaxes(type="log", row=1, col=1)
# plot1_2

8.5 Unbalance response

ROSS’ method to simulate the reponse to an unbalance is run_unbalance_response(). This method returns the unbalanced response in the frequency domain for a given magnitide and phase of the unbalance, the node where it’s applied and a frequency range.

ROSS takes the magnitude and phase and converts to a complex force array applied to the given node:

\[\begin{split}force = \left(\begin{array}{cc} F \cdot e^{j\delta}\\ -jF \cdot e^{j\delta}\\ 0\\ 0 \end{array}\right)\end{split}\]

where: - \(F\) is the unbalance magnitude; - \(\delta\) is the unbalance phase; - \(j\) is the complex number notation;

When plotting the results, you can choose to plot the: - Bode plot options for a single degree of freedom: - amplitude vs frequency: plot_magnitude() - phase vs frequency: plot_phase() - polar plot of amplitude vs phase: plot_polar_bode() - all: plot() - Deflected shape plot options: - deflected shape 2d: plot_deflected_shape_2d() - deflected shape 3d: plot_deflected_shape_3d() - bending moment: plot_bending_moment() - all: plot_deflected_shape()

run_unbalance_response() is also able to work with clustering points ( see section 7.4.1 ).

8.5.1 Running unbalance response

To run the Unbalance Response, use the command .unbalance_response()

In this following example, we can obtain the response for a given unbalance and its respective phase in a selected node. Notice that it’s possible to add multiple unbalances instantiating node, magnitude and phase as lists.

The method returns the force response array (complex values), the displacement magnitude (absolute value of the forced response) and the phase of the forced response.

Let’s run an example with 2 unbalances in phase, trying to excite the first and the third natural vibration mode.

Unbalance1: node = 29
            magnitude = 0.03
            phase = 0
Unbalance2: node = 33
            magnitude = 0.02
            phase = 0

n1 = 29
m1 = 0.003
p1 = 0

n2 = 33
m2 = 0.002
p2 = 0

frequency_range=np.linspace(315, 1150, 101)
results2 = rotor3.run_unbalance_response([n1, n2], [m1, m2], [p1, p2], frequency_range)

8.5.2 Plotting results - Bode plot

To display the bode plot, use the command .plot(probe)

Where probe is a list of tuples that allows you to choose not only the node where to observe the response, but also the orientation.

Probe orientation equals 0º refers to +X direction (DoFX), and probe orientation equals 90º (or \(\frac{\pi}{2} rad\)) refers to +Y direction (DoFY).

You can insert multiple probes at once.

# probe = (probe_node, probe_orientation)
probe1 = (15, 45) # node 15, orientation 45º
probe2 = (35, 45) # node 35, orientation 45º

results2.plot(probe=[probe1, probe2], probe_units="degrees")

# converting the first plot yaxis to log scale
# plot2 = results2.plot(probe=[probe1, probe2], probe_units="rad")
# plot2.update_yaxes(type="log", row=1, col=1)
# plot2

8.5.3 Plotting results - Deflected shape

To display the deflected shape configuration, use the command .plot_deflected_shape()

results2.plot_deflected_shape(speed=649)

8.6 Time response

ROSS’ method to calculate displacements due a force in time domain is run_time_response(). This function will take a rotor object and plot its time response given a force and a time array.

The force input must be a matrix \(M \times N\), where: - \(M\) is the size of the time array; - \(N\) is the rotor’s number of DoFs (you can access this value via attribute .ndof).

Each row from the matrix represents a node, and each column represents a time step.

Time Response allows you to plot the response for: - a list of probes - an orbit for a given node (2d plot) - all the nodes orbits (3d plot)

8.6.1 Running time response

To run the Time Response, use the command .run_time_response().

Building the force matrix is not trivial. We recommend creating a matrix of zeros using numpy.zeros() and then, adding terms to the matrix.

In this examples, let’s create an harmonic force on node 26, in \(x\) and \(y\) directions (remember index notation from Frequency Response (section 7.4.2). We’ll plot results from 0 to 10 seconds of simulation.

speed = 600.0
time_samples = 1001
node = 26
t = np.linspace(0, 16, time_samples)

F = np.zeros((time_samples, rotor3.ndof))

# component on direction x
F[:, 4 * node + 0] = 10 * np.cos(2 * t)
# component on direction y
F[:, 4 * node + 1] = 10 * np.sin(2 * t)

response3 = rotor3.run_time_response(speed, F, t)

8.6.2 Plotting results

There 3 (three) different options to plot the time response: - .plot_1d(): plot time response for given probes. - .plot_2d(): plot orbit of a selected node of a rotor system. - .plot_3d(): plot orbits for each node on the rotor system in a 3D view.

Ploting time response for list of probes

probe1 = (3, 0)   # node 3, orientation 0° (X dir.)
probe2 = (3, 90)  # node 3, orientation 90°(Y dir.)

response3.plot_1d(probe=[probe1, probe2], probe_units="degree")

Ploting orbit response for a single node

node = 26
response3.plot_2d(node=node)

Ploting orbit response for all nodes

response3.plot_3d()

8.7 Undamped Critical Speed Map (UCS)

This method will plot the undamped critical speed map for a given range of stiffness values. If the range is not provided, the bearing stiffness at rated speed will be used to create a range.

Whether a synchronous analysis is desired, the method selects only the foward modes and the frequency of the first forward mode will be equal to the speed.

To run the UCS Map, use the command .plot_ucs().

8.7.1 Running and plotting UCS Map

In this example the UCS Map is calculated for a stiffness range from 10E6 to 10E11 N/m. The other options are left to default.

stiff_range = (6, 11)
ucs_fig = rotor3.plot_ucs(stiffness_range=stiff_range, num=20, num_modes=16)
ucs_fig